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I Abstract: We study the /i < branch of the minimal supergravity ansatz of the 

Q^l minimal supersymmetric standard model. The extent to which /i < is disfavoured 

p i' compared to ;U > in global fits is calculated with Markov Chain Monte Carlo 



X 



methods and bridge sampling. The fits include state-of-the-art two-loop MSSM con- 
tributions to the electroweak observables Mw and sin^ 6'^, as well as the anomalous 
magnetic moment of the muon {g — 2)^, the relic density of dark matter and other 



\-{ ■ relevant indirect observables. /i < is only marginally disfavoured in global fits and 

should be considered in mSUGRA analyses. We estimate that the ratio of probabil- 
ities is P(/i < 0)/P(/i > 0) = 0.07 - 0.16. 
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1. Introduction 

There has been increasing recent attention on global fits of various indirect data 
to minimal supergravity (mSUGRA), also sometimes called the constrained min- 
imal supersymmetric standard model (CMSSM) |[l|, [^, |[, ^. mSUGRA makes 
phenomenological analysis of the minimal supersymmetric standard model (MSSM) 
tractable via the low number of free parameters. In fact, the scalar masses mo, gaug- 
ino masses M1/2 and trilinear coupling Aq are assumed to be universal at a gauge 
unification scale Mqut ~ 2 x 10^^ GeV. If the MSSM is present in nature and if 
the mSUGRA universality assumptions are approximately correct, chi-squared or 
probability distributions for potential collider /dark matter observables can be de- 
rived. Early fits [I|, [^, ^ Q necessarily had fixed input parameters to reduce the 
dimensionality of the input parameter space, making scans practicable. It is usually 
assumed that neutralinos constitute the current cold dark matter content of the uni- 
verse, since they are weakly interacting, electrically and colour neutral and stable. 
The predicted value of dark matter relic density floMh'^ is a very strong constraint 
on viable mSUGRA parameter space, effectively reducing its dimensionality by 1. 

The accuracy of the inferred value of QoMh"^ from WMAP data makes a global fit 
to all of the relevant mSUGRA parameters potentially difficult because the system is 
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rather under-constrained, possessing narrow, steep valleys of degenerate minima. 
If the MSSM is confirmed in colliders, it will hopefully be possible to break such 
degeneracies with collider observables. This does not help us at present, where we 
want to provide a sort of 'weather forecast' for future colliders. It was indicated 
in ref. that the powerful Markov Chain Monte Carlo (MCMC) technique might 
allow us to find the probability distribution of a fully global fit to indirect data. 
Two of us went on |^ to demonstrate that MCMCs do indeed allow such a fit, 
investigating collider observables. One of us examined the effect of a naturalness 
prior [§]. Our results were confirmed and expanded in Ref. 0, also utilising the 
MCMC method and including a one-loop MSSM calculation of the W-boson mass 
M\Y and the weak leptonic mixing angle sin^ 6'^ in the likelihood density. The purpose 
of MCMC mSUGRA global fits is two-fold: as well as producing interesting and 
useful physics results in themselves, we may profit from the experience of utilising 
and developing the MCMC tools, which could prove very useful when analysing 
future collider data. 

It is our purpose in the present paper to extend the previous /i > global fits 
to < 0. Besides the observables studied in we now also include Mw and sin6'^ 
in our analysis, as done in e.g. ^. Being highly sensitive to new physics these 
very accurately measured quantities play a key role in the electroweak sector and are 
therefore also of great interest when it comes to further constraining the mSUGRA 
parameter space. It was shown in the literature and that the one-loop predictions for 
the two observables alone do not bear enough accuracy to make reliable predictions. 
In fact, the pure one- loop predictions can lead to results contradictory to the state- 
of-the-art predictions |]IU| used in our analysis. These contain the known higher 



order contributions from both the Standard Model and the MSSM. Extending our 
analysis to negative values of /i it is crucial to further use a very accurate prediction 
for [g — 2)^. The dominant two-loop corrections [jlll to this quantity are therefore 
also taken into account in the present analysis. It is well known that the measured 
anomalous value of the magnetic moment of the muon [g — 2)^ is roughly 2a above 
the Standard Model (SM) predicted value. This positive contribution is predicted 
by some regions of mSUGRA parameter space, provided /x > 0. The "dark side" of 
mSUGRA (i.e. /i < 0) provides a negative contribution, thereby being disfavoured by 
the {g — 2)^ measurement. In a global fit, one can trade likelihood penalties between 
different observables and the conclusion that /i < is disfavoured to roughly 2cr is 
not at all obvious. We will calculate the extent to which the dark side is ruled out by 



using MCMCs with "bridge sampling" |]12|. We will encounter problems associated 



with isolated likelihood density maxima in the dark side, potentially ruining MCMC 
convergence. Fortunately, bridge sampling provides a solution to the convergence 
issue and we are able to calculate the degree to which the dark side is disfavoured 
with respect to yU > 0. As well as extending previous analyses to < 0, we have 
made several technical improvements in the calculation of the likelihood compared 
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with previous attempts in the hterature. 

If the hghtest supersymmetric particle decays into Standard Model particles, as 
is the case in R-parity violation for instance, its relic density will be essentially zero 
today. In that case, one requires to obtain the WMAP fitted flDMh"^ from some 
other source than neutralinos (gravitinos or hidden sector matter for instance). In 
order to investigate this case, we will also perform the fits for the case where all 
relevant data except QoAih"^ are included in the likelihood density. Such fits will 
help us to understand the impact of VLoKih'^ in constraining the model, as well as 
being relevant for the R-parity violating mSUGRA in the limit of small R-parity 
violating couplings. 

We now go on to detail the various constraints used on the model in section |[ 
The results of the dark side fitting procedure are compared and contrasted against the 
better known /i > ones in section before the effect of dropping the dark-matter 
constraint is examined in section |[ Closing remarks are presented in section |^. A 
presentation of the fitting procedures is confined to the appendix: Markov Chain 
Monte Carlos and bridge sampling are discussed in appendix |^. Convergence prob- 
lems and their resolution are discussed in appendix 



2. Constraints 



We vary 8 input parameters rele- 
vant to the model. The range of intrin- 
sically mSUGRA parameters consid- 
ered is shown in Table |l], where tan (3 
is the ratio of the two MSSM Higgs 
doublet vacuum expectation values. We 



mSUGRA parameter 


range 




-4 TeV to 4 TeV 


mo 


60 GeV to 4 TeV 


Ml/2 


60 GeV to 2 TeV 


tan (5 


2 to 62 


SM parameter 


constraint 




127.918±0.018 




0.1176±0.002 




4.24±0.11 GeV 


nit 


171.4±2.1 



use Ref. for the QED coupling con- 
stant a*^"^, the strong coupling con- 
stant af^{Mz) and the runnmg mass 
of the bottom quark mf,(mfe)*^'^, all in 
the MS renormalisation scheme The 
recent Tevatron top mass mt measure- Table 1: Input parameters 

ment |]T5[ is also employed. These SM 

inputs are shown in Table |I[ Experimental errors are so small on the mass of the 
Z° boson Mz and the muon decay constant that we fix them to their central 
values of 91.1876 GeV and 1.16637 x 10"^ GeV'^ respectively. The Standard Model 
(SM) input parameters are allowed to vary within 4(T of their central values but a 
penalty 

{Ci-pi{m)f 



'X.i 



(2.1) 



3 



is applied for observable i. Ci denotes the central value of the experimental measure- 
ment, pj(m) represents the value "predicted" at any stage of the MCMC sampling 
given knowledge of the model m presumed to be "true" at that point. Finally is 
the standard error of the measurement. Equivalently, expressing this in the language 
of likelihoods, we are assuming that each of these measurements have Gaussian er- 
rors,^ and that the likelihood distribution Ci = p(cj|m) for any one measurement 
may be written in the following way: 

Ci = p(ci|m) = 



27r cr. 



:exp [-X.V2]. 



[2.2) 



The normalisation constant \l 2-1:01 may be ignored in subsequent calculations as 
the absolute value of Ci will never be needed. It will only be necessary to know the 
ratios of values of Ci at neighbouring points in the MCMC chain, or between chains 
in which the neglected constants are identical. 

In order to calculate predictions for observables from the inputs in Table |T], we 



use S0FTSUSY2 . . 7 fTl to first calculate the MSSM spectrum. We apply the bounds 





37 




67.7 


nig 


195 


ruf^ 


76 


mr 


88 




86.4 




91 


^gR 


250 




43.1 















Table 2: Lower bounds applied to sparticle mass predictions (in GeV). 

in Table ^ in order to take into account 95% limits coming from negative sparticle 
searches fl^. Any point transgressing these bounds is given a zero likelihood density 
(or, equivalently, an infinite x^)- Also, we set a zero likelihood for any inconsistent 
point, e.g. one which does not break electroweak symmetry correctly, or a point that 
contains tachyonic sparticles. For points that are not ruled out, we then link the 
MSSM spectrum via the SUSY Les Houches Accord ^ to micrOMEGAsl . 3 . 6 [|1|], 
which then calculates fljjMh'^, the branching ratios BR{b 57) and BR{Bs — > 
and {g - 2)^. 

The measured value of the anomalous magnetic moment {g — 2)^ is in confiict 
with the SM predicted value by 



(5a„ = 6 



(9-2) 



(22± 10) X 10^ 



10 



:2.3) 



This excess may be explained by a supersymmetric contribution, the sign of which 
is identical to the sign of the superpotential /i parameter ||19|. After obtaining the 



^The experimental constraints on BR{Bs and the LEP constraints on the Higgs mass, 

each described later, are not Gaussian constraints and must therefore be treated differently. Nothing 
prevents us from continuing to parametrise their likelihood distributions in the same way, however, 
but it should be realised that a consequence of this is that the "x-squared penalty" (i.e. — 21og£i) 
will not be parabolic, as Ci is not a Gaussian distribution in these cases. 
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one-loop MSSM value of {g — 2)^ from micrOMEGAs, we add the following dominant 
2-loop corrections |Tl|, the logarithmic piece of the 2-loop QED contribution, 



16 



two- loop stop-higgs and chargino-stop/sbottom contributions. 

The Tevatron has recently been 
instrumental in bounding the branch- 
ing ratio of the rare decay channel 14 
Bs fJ'^fJ'^ mi- Such bounds help 



constrain the mSUGRA parameter space 
We apply a penalty on the value 
predicted by micrOMEGAsl . 3 . 6 de- 
rived from CDF Tevatron Run II data 
The resulting penalty is shown in Fig. |I|. 

Recently, it has been claimed that 
light sparticles are preferred by the 
two weak observables sin^ 6'^ and Mw 
[|, In Ref. [|, Mw and sin^^^ 




-7.4 -7.2 -7 -6.8 
logio [BR(B,->nV)] 



6.4 



from 



were used at one loop order to help 
constrain mSUGRA in a global MCMC 

fit. The preference for such hght SUSY Figure 1: penalty on BR(5, 
was not particularly evident in the 
global fits. We examine the mSUGRA 

predictions for Mw and sin^6'^(eff) in Figs. | and ^ for Aq = 0, tan/3 = 10, yU > 0, 
equal rriQ and M1/2 and central experimental values for the other inputs. For the 
"SOFTSUSY" lines, the default SOFTSUSY calculation is used. This contains the full 
SOFTSUSY MSSM contributions to the leptonic mixing angle sin^^^ and Mw- It 
also contains the dominant 2-loop Standard Model contributions to Mw- For the 
lines marked "2-loop" , the SUSY Les Houches Accord is used to communicate with a 
currently private code that calculates the W-boson mass Mw and the effective 
leptonic mixing angle variable sin^ 6*1^, calculated to two loops in the dominant MSSM 
parameters. We use the most general MSSM result for the full one-loop contribu- 
tions. Besides all known corrections due to SUSY particles, the full SM contributions 
are also included in the predictions for Mw and sin^ 6^^, leading to the currently most 
accurate predictions within the MSSM. The "SM" lines show the SM limit, where 
all corrections involving sparticles are dropped, i.e. the state-of-the-art SM results 
p5| , with MffSM = Mh- They vary slightly with mo = M1/2 because the varying 
mSUGRA parameters produce different values of the Higgs boson mass nih. The 



horizontal lines on the figures show the current la experimental limits P7| , 



w 



80.392 ± 0.031 GeV, 



sin^ el = 0.23153 ± 0.00020, 



(2.4) 



where we have added experimental and theoretical errors in quadrature. 
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SOFTSUSY 
2-loop 
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expt 1 sigma 
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0.232 



0.2318 



crO.2316 
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exp 1 sigma 
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mo=Mi/2 (TeV) 



Figure 2: Various approximations to Figure 3: Various approximations to 



Mw in mSUGRA 



sin^ e,n in mSUGRA. 



The theoretical errors in the predicted Myy and sin^ 6'^ are estimated to be 10 
MeV and 12x10"^ respectively 0. We use these uncertainties for the purposes of 
comparison, although they have been slightly reduced recently by the addition of 
additional two-loop corrections taken into account in the present analysis |10|, ^ . 
We see from the SOFTSUSY line in the figure that the prediction of does not 
have a strong preference for light SUSY, since the model is within the la errors 
up until niQ = M1/2 = 4 TeV. In actual fact, only very light SUSY masses are 
disfavoured by the "SOFTSUSY" line, leading to predictions above the la-range. 
The situation is similar for the sin^ 0^ "SOFTSUSY" line, where again only very 
light SUSY masses lead to predictions outside the la-range. However, using the best 
available predictions, corresponding to the "2-loop" lines, a preference for light mo = 
Ml/2 can be seen in the prediction for M^. The SM curve which lies just below the 
la- interval is approached from above in the decoupling limit, furthermore indicating 
a slight preference of the MSSM over the SM. The preference for light SUSY is 
not as striking for sin^^^. Here most of the ttiq = M1/2 values are doing equally 
well, which is mainly due to the fact that the SM prediction for sin^ 6'^ is already 
well within the la-range. With the behaviour of the "one-loop" curve and the best 
available result being qualitatively different, it is desirable to use the more accurate 
result for Mw and sin^ 6*^ when calculating the contributions of Mw and sin^ 
using Eq.p.l|. Although Ref. only used the one-loop predictions, the theoretical 
errors were correspondingly enlarged in order to take the larger uncertainty from 
higher order terms into account. 
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LEP2 constraints on the lightest CP- 
even hieffs mass are included further 
likelihood penalty following a parameteri- 



sation of LEP2 data in the SM limit 30 



For the LEP2 constraints, the SM limit 
is a good approximation for mSUGRA, 
since sparticle mass limits imply that we 
must be near the decoupling regime of 
the MSSM 0. We estimate that the 



S0FTSUSY2 .0.7 determination of mt has 
a 2 GeV theoretical error in mSUGRA, al- 
though it may be somewhat larger in the 
general MSSM [^. We therefore smear 
the parameterised LEP2 Higgs likelihood 
density CLEP2{fnh) with a Gaussian dis- 
tribution of width CTft = 2 GeV: 




110 112 114 116 118 120 122 124 
mj^o (GeV) 

Figure 4: LEP2 higgs penalty paid 



dx 



2TTah 



(2.5) 



The result of this procedure leads to the effective Ax^ = —2 In Ch penalty shown in 
Fig. ^ The slight excess of candidate Higgs events over the background prediction 
at LEP2 can be seen by a negative Ax^ penalty in the figure for rrih ~ 116 — 121 
GeV. 



The rare bottom quark branching ratio is BR{b 57) is constrained to be ^3 



BR{h 57) = (3.55 ± 0.38) x 10" 



(2.6) 



obtained by adding the experimental error with the estimated theory error of 
0.3x 10~^ in quadrature. Very recent estimates ^ of BR{h — > 57) are compatible 



BR{b 37) is substituted for pi in Eq. ^]T] in order to calculate 



with Eq. |2.6| at the Icr level, although the error has decreased. Our prediction of 



(2.7) 



We use the WMAP3 fS^] power law A-CDM fitted value of 

= ilDMh'' = 0.104+°-°°^3 



-0.0128 



for the dark matter relic density of the universe. We initially assume that the neu- 
tralinos are stable and that they constitute the whole of the dark matter relic density. 
Eq. ^A] is used to calculate Xq; with = 0.0073 for a prediction lower than the 
central experimental value and 0.0128 otherwise. 

Having described the calculation of the likelihood associated with each individual 
measurement, we are now in a position to define the likelihood of the set of all 
measurements or observables, taken together. We are required to calculate the joint 
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(total) likelihood C of all the measurements given the truth, i.e. in the notation of 
Eq |2]l| we want to know: 



L = p(measurements|true model) (2.8) 

= p(ci,C2,...|m) (2.9) 

= p(ci|m) ■ p(c2|ci, m) ■ p(c3|ci, Ca, m) ■ . . . (2.10) 

= £i ■ p{c2\ci, m) ■ p{c3\ci, C2, m) ■ . . . (2.11) 

which is not in general equal to 

Ci-C2-C3-... (2.12) 

unless we can be confident that for each measurement 

p(c3|ci,C2,m) ^p(c3|m) etc. (2.13) 

Fortunately we can be confident that Eq |2.13| does hold in our situation, as we have 
deliberately constructed m to be broad enough such that all fundamental parameters 
which might reasonably be expected to correlate any two of the measurements are 
included within it^. We are therefore at liberty to write: 

£ = £1 ■ £2 • A ■ • • • (2.14) 

= e-^»^?/2 (2.15) 

for the total likelihood, and we can be confident that this product takes into account 
the expected correlations between all the observables contained, due to the nature 
of the space {m} of models considered. 



3. Dark Side Fits 

We now compare and contrast the dark side fits to those with n > 0. In Figs, ^(a) 
and|^(b), we show the posterior probabilities P marginalised^ to the mo — M1/2 plane 
for both signs of fi. As with all 2-dimensional marginalised plots in this paper, we 
bin the plane into 75x75 2-dimensional bins. The colour bar on the right hand side 
of the figures shows the posterior probability P of each bin divided by the maximum 



^If the design of m were not broad enough, then m would have to be extended. For example: 
were it the case that the up quark mass niu was expected to significantly correlate two or more of 



the observables, then for Eq 2.13 to continue to hold, rn„ would have to be added to m enlarging 
its dimension by one. 

■^For readers unfamiliar with the term: marginalisation means "integrated over the unseen di- 
mensions of parameter space" in this context. 
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Figure 5: mSUGRA Fits for (a) ;U < (b) /j, > marginalised to the mQ-Mi/2 plane. 
The posterior probability is indicated by the bar on the right hand side. The inner (outer) 
contours show the 68% (95%) confidence region respectively. 



posterior probability of any bin in the plot. In Fig. ^(a), the only 68% contour^ is 
at the lowest values of mo: all other contours are 95% confidence region contours 
due to the lowish likelihood densities. The brokenness of the contours is a result 
of statistical fluctuations in the results. Although these are visible, there is clearly 
reliable information their trajectories in the plot. The n < plot displays two isolated 
local maxima, whereas the 95% confidence region of yU > is continuously connected. 
We discuss the relative normalisation of the two /i < maxima in appendix |B[ 
In the 95% confidence region of Fig. |^a closest to the origin, the relic density is 
dominantly depleted by either stop co-annihilation ^ ^ ixi tg or stau 
co-annihilation TiXi ~^ '^1- the other hand, the 95% region at higher rriQ — 



Ml/2 consists of resonant Higgs annihilation regions P2| , ^ where XiXi ~^ 
h, — > bb/T~^T~ and the focus point region where the LSP has a significant higgsino 
component and XiXi ~^ WW, tt P5| , 
annihilation 



I, ff^. There was no significant stop co- 



region for /i > 0. 

We now include some other mar ginalisat ions on the other 2-dimensional parame- 
ter planes for yU < mSUGRA with a fiat prior. They are displayed in Figs. P(a)-(d). 
Figs. |](a) and ^(b) show that the probability density for the 2-chain co-annihilation 
sample is not separated from the other sample in either the M1/2 — Aq plane or the 
tan P-Ao plane (the almost disconnected region at the bottom of Fig. |](a) consists of 



■^Note that the confidence regions in Figs. || (and those in later plots) should strictly be referred to 
as "Bayesian credible intervals" (each region contains a fixed amount of the posterior probability) to 
distinguish them from the related concept from Frequentist Statistics called a "confidence interval" . 
Usage of the term "credible interval" is not common in High Energy Physics, however, and we will 
stick to "confidence region" . 
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H 



(c) 



> 
H 




- 0.' 



Figure 6: Constraints from global fits with /i < mSUGRA marginalised to 2-dimensional 
parameter planes. We have assumed a flat prior. The posterior probability is indicated by 
the bar on the right hand side. The inner (outer) contours show the boundary of a 68% 
(95%) confidence region. 



the light h^-pole region). There is only a modest separation in the Mi/2-tan/5 plane: 
the 68% contours do not connect the two regions, whereas the 95% contours do. 
Fig. |](d) shows that the connection between the two samples in the mo — tan (3 plane 
is marginal. The rriQ — Aq marginalisation was useful for investigating the physics 
behind the two isolated probability maxima, and is displayed in appendix p. 

The probability distributions of the masses of selected MSSM particles are shown 
in Fig. 1^ for both signs of /i in mSUGRA. The /z < sample has not been normalised 
to the correct relative normalisation compared to yii > in the figure. The lightest 
CP-even higgs, the gluino and the lightest neutralino all have mass distributions 
that are remarkably similar for either sign of /i. However, in Fig. 0(c), we see that 
the fi < sample has a flat plateau for higher squark masses, whereas the yU > 
sample tails off somewhat. High squark masses will result in smaller total SUSY 
cross-sections at the LHC but the lighter gluino should still provide enough events 
for discovery if rrig < 2 TeV [^, Sharp peaks at low gluino and neutralino 
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Figure 7: Probability distributions in mass of (a) the lightest CP even higgs, (b) gluino, 
(c) the left-handed squark and (d) the neutralino. Flat priors have been assumed. 

masses are due to the /i°-resonance annihilation region 0. The broader peak of the 
/i < curve in Fig. ^(c) is due mostly to the co-annihilation sample. The significant 
probability densities for large gluino and squark masses are rather alarming, as then 
SUSY detection at the LHC would require a longer running time. Gauginos are 
not so sensitive to the range of the prior in Table ^ but the sfermions are and a 
reduced range makes lighter sfermions more likely. Also, naturalness priors have 
a large impact, reducing the likely sparticle masses. 

3.1 Further investigations of the fits 

We examine the best-fit points from each sampling in Table 0. We see from the table 
that, as expected, the fi > sample has a better best-fit point and a correspondingly 
lower x^, mainly due to the much better fit to 5a^. In agreement with Refs. [^, 0], 
the best-fit /i > point is at light SUSY masses. The mo and M1/2 parameters are 
smaller for the /i > case than for /i < 0, corresponding to lighter sparticles and 
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fi<0 


/i > 




< 


> 


mo/GeV 


3610 


156 


^a^/lO-io 


-0.4 


14.2 


Mi/2/GeV 


93 


569 


BR{b^ s-f)/ 10'^ 


3.65 


3.41 


Ao/GeV 


-56 


270 


BR{Bs /i+/i~)/10~9 


3.1 


3.7 


tan f3 


6.0 


24.1 


sin^ ^^(eff) 


0.23153 


0.23152 




0.102 


0.101 


Mw/GeV 


80.382 


80.368 


m/i/GeV 


117.5 


115.8 




4.5 


1.5 



Table 3: Best-fit points from the MCMC samplings for each sign of /i. = Yli xf 
the Standard Model inputs are close to their experimental central values in each case. 



the larger contribution to the anomalous magnetic moment of the muon. If we take 
the n < best-fit point and flip the sign of fi, we find that the point does not break 
electroweak symmetry properly and is excluded. The fi > best-fit point has a total 
of 181 when the sign of fi is flipped, mostly due to an increase in the predicted 
value of QDAih"^ to 0.195. 

The region of smaller mo, M1/2 is more probable for /i > than for /i < 0: this 
will lead to a relatively heavier /i < spectrum. Our fi > results are generally 
similar to previous analyses which did not include M^y and sin^ 6'^ as constraints in 
Ref. and to those which included the one-loop SOFTSUSY prediction for Myy, sin^ 6*^ 
with enlarged theoretical errors 0. This seems in contradiction to the conclusions 



of Ellis et al H, 24|, where it is claimed that the electroweak variables prefer light 



SUSY. Indeed, Fig. ^ indicates that My/, sin^ 6*^ do mildly prefer light SUSY but 
our results show that this preference is washed out in the global fits. Our results 
allow for heavier sparticles than Ellis et al, mainly because we have chosen to allow 
more relevant parameters to vary: 8 compared to 2 in their paper (one dimension is 
fixed by requiring the relic density prediction to be the central WMAP-constrained 
value). Their fits are for different discrete values of fixed tan/5, but if it were allowed 
to vary, we believe that the confidence regions there would be enlarged. In the 
present paper, we also obtain additional smearing from allowing m^, a^, a and mj, to 
vary. Figs. ||(a) and |(b) show the probability distributions for My/, sir? 6*^ coming 
from the mSUGRA fits for both signs of /i, although the sign does not make much 
difference. We see that the My/ prediction coming from the fits is skewed towards 
values lower than the central empirical value, corresponding to a preference for heavy 
SUSY from the rest of the fits. Fig. ^ confirms that heavier SUSY tends to have 
lighter values of My/- From Fig. ^ we see that heavy SUSY tends to be on the upper 
la empirical limit of sin^^^. Fig. H(b) does show evidence for this skew, which is 
rather mild. 

The strongest constraint in the fits comes from the dark matter relic density. In 
Fig. 1^, we show the probability densities resulting from the fits for the dark matter 
relic density. Each curve is normalised slightly differently to allow better viewing of 
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Figure 8: Posterior probability distributions for weak observables in mSUGRA. 
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the results. We see that both the /i > and the < VtoMh^ distributions follow 
the empirical constraint closely, except for a slight excess just lower than the central 
value. 

Another important aspect influencing the fits is the integrated volume of the 
probability density, which is automatically taken into account in a Bayesian analysis, 
as was demonstrated and pointed out in Ref . @ . The usual arguments based purely 
on values of above the best-fit value are valid when the probability density function 
is Gaussian in the interesting parameters (which is certainly not the case here, as 
even a cursory glance at Fig. ^ allows). Thus even if, say, the stop co-annihilation 
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|i<0 



region had a much lower than other regions of the fits,^ the fact that its volume 
in 8-dimensional input parameter space is much smaller than the other regions will 
automatically disfavour it since its integrated probability will be low. 

The LEP2 higgs constraints shown in 
Fig. ^, along with the current empirical 
value of rrit shown in Table |T] favours rather 
heavy mSUGRA. Indeed, Fig. ITUI shows 
that rrit is skewed to somewhat higher masses 
than the empirical constraint, illustrating 
this tension (again, we have altered the 
normalisation of the curves slightly for clar- 
ity). Inspection of the a{Mz)^, a,{MzY^, 
mb{mb) inputs show that they follow their 
empirical probability distributions very closely. 
There is a large volume of parameter space 
for the dark matter annihilation region 

at high tan/? > 10 particularly for u > 0, t-,. o u u-r+ a- + -u +• t 

° ^ ^ J r- ) l< igure 11: rrobability distribution tor 

as shown in Fig. 0. High values of tan/? ^^^^ mSUGRA 
mean that light SUSY is disfavoured by 

BR[b — > 57] and BR[Bs — > since the data disfavours SUSY contributions 
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which are approximately proportional to ta.ia^ P/Mg^gy and tan^ /?/M|^^y respec- 
tively. The distributions of these two observables are shown in Fig. |T2|. As can be 
seen from the figure, the sign of the SUSY contribution to each observable depends 
upon the sign of /i. The maxima in each figure correspond to observables being close 
to their SM limit. BR[b 57] prefers /i > mildly, as /i < tends to predict a 
BR[b — i> 57] larger than the central empirical value. 

is expected to be the observable that most strongly discriminates between 
the two signs of /x. We plot its distribution in each case in Fig. |l3|(a). Since < 
has the wrong sign of compared to experiment, the probability density bunches 
around zero. Clearly heavy SUSY with a less negative SUSY contribution is favoured. 
However, we see that the /i > distribution also prefers smaller SUSY contributions 
than the data in the global fits. This was initially unexpected, and will lead to 
/i < mSUGRA being less ruled out. In order to understand this behaviour better, 
we re-weight the /i > sample in various ways. Sa^ is approximately proportional 
to taii^ P/Mg^gY need large tan/? and rather light SUSY in order to get a 

sizable value in line with the central experimental value. As explained above, this is 
somewhat in conflict with LEP2 Higgs constraints and BR{b 57). We re- weight 
the fi > chains, dividing by a number that removes the likelihood contribution 
from the LEP2 Higgs constraint and the BR{b 57) measurement, i.e. ChCBR{b^s'y)- 
The probability distribution of 6a^ resulting from this procedure is marked as the 
"reduced" curve in Fig. p^ (b). It e xtends to somewhat higher and more central values 
^In reality, the stop co-annihilation region fits the data rather marginally. 
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Figure 12: Probability distributions of (a) BR{b — > 57) and {h)BR{Bs — > //^/U ) in 
mSUGRA. 
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Figure 13: Probability distributions for in mSUGRA (a) for flat priors and both signs 
of fi and (b) for ^ > and various different priors. The key is explained in the text. 



of 5a^, showing that some of the skew in the (5a^ distribution does come from the 
LEP2 Higgs and BR{b —>■ s'j) measurements. However, the effect is rather small and 
the resulting distribution is rather far from the experimental distribution, indicating 
a further effect. There could be a significant volume effect if regions of parameter 
space that have central values of 6a ^ have a small volume. This would render our 
results sensitive to the prior, since by changing the measure of the input parameters, 
we can change the volume measure 0. In order to investigate the effects of this, we 
re-weight the fi > chains by a factor l/(moMi/2)- Such a re- weighting mimics the 
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effect of using a logarithmic prior on mo and M1/2 since 

/f , , (ilnmn (ilnMi/2 
P{mo,Mi/2) rflnmo rflnMi/2 = / P(mo,Mi/2)-^ ^ c/mo c?Mi/2 

y ctmo ctiWi/2 

One might consider such a prior on technical naturalness grounds. The "log prior" 
results are displayed in Fig. |TB|(b). They have a fatter tail out to higher values of 
Sttfj, than the fiat prior sample or the "reduced" sample. When we perform the re- 
weighting in Eq. as well as the one to remove BR{b 57) and LEP2 Higgs 
constraints, we obtain the "log prior, reduced" curve in the figure. This shows a yet 
fatter tail and is starting to approach the empirical probability distribution imposed 
on the fits. Our results are obviously somewhat dependent upon the prior. This is 
essentially because there is not enough precise data yet to constrain mSUGRA very 
strongly. We must bear in mind this dependence upon the prior, and investigate 
different priors when we estimate P_/P+ below. 

The total normalisation of the fi > 0, ^ < 
samples is shown in Table ^ for various different pri- 
ors. The "small" prior is a fiat prior with a reduced 
range compared to the one displayed in Table 0. We Table 4: Ratios of integrated 
filter the chains to only include points with mo < 2 Probability for different signs of 
TeV and |y4o| < 2 TeV. The table shows that /i < ^• 
is somewhat disfavoured for the smaller ranges and 

more disfavoured still for the log prior. The result is somewhat sensitive to the 
prior, indicating the need for more data. We therefore prefer to quote a range for 
P_/P+ = 0.07 — 0.16 depending upon the prior. This range is a focal result of the 
present paper. 



prior 


flat small log 


P-/P+ 


0.16 0.12 0.07 



4. Fits Without WMAP3 

We now briefly examine the effect of removing the dark matter penalty from the 
fits. We could in principle re-weight the chains in order to do this, but we find 
that that leads to statistical fluctuations in the results that are too large. Initial 
investigations revealed that the efficiency becomes much higher when we remove 
the dark matter relic density contribution to the total x^- We are able to increase 
the widths of the proposal distribution to 197 GeV for mo, 100 GeV for M1/2, 400 
GeV for Aq, 8.5 for tan/? and la for the Standard Model inputs and still achieve 
an efficiency of around 35%. This has the consequence that the chains explore the 
parameter space much quicker than in the previous section, and so less MCMC steps 
are needed. We run a further 9x200 000 MCMC steps for each sign of fi. This time. 
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Figure 14: Constraints from global fits without the dark matter constraint with (a) ;U > 
and (h) fj, < mSUGRA marginalised to the mQ-Mi/2 plane. We have assumed a flat prior. 
The posterior probability is indicated by the bar on the right hand side. The inner (outer) 
contours show the boundary of a 68% (95%) confidence region. 
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both signs of /i have excellent convergence statistics, R being different to 1 at the 
per-mille level only in each case. 

We see from Fig. that removing the 
dark matter constraint yields a very different 
picture for the mSUGRA probability distri- 
butions in the mo-Mi/2 plane. This confirms 
our statement that many of the features seen 
in the previous section are due precisely to 
that constraint. The probability distribution 
is now much flatter in the input parameter 
space. The disallowed region at low mo and 
high Ml/2 is due to the no-charged LSP con- 
straint. The disallowed region at low rriQ and 
low Ml/2 is due mostly to the combined ef- 
fect of the {g — 2)^, BR{b 57) and rrih 
constraints. The disallowed region is signifi- 
cantly larger for /i < than for /i > 0, due 

to those three observables. Marginalisations in other mSUGRA input parameter 
planes tell a similar story, more featureless than the fits including the dark matter 
constraint. As mentioned in the introduction. Fig. |14| covers the case of R-parity 
violating mSUGRA when the R-parity violating couplings are smaller than about 
0.1. For larger R-parity violating couplings, one would have to include them in the 
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Figure 16: Jl/i^ distributions without 
imposing the dark matter constraint. 
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Figure 15: Probability distributions in mass of (a) the lightest CP even higgs, (b) gluino, 
(c) the left-handed squark and (d) the neutralino. Flat priors have been assumed. 



renormalisation group equations to obtain accurate results. 

The more featureless fits have predictable effects on the mass spectrum of the 
MSSM, as shown in Fig. |l^. We have fixed the ranges of the abscissas to be identical 
to those in Fig. ^ in order to facilitate comparison. The lightest CP-even Higgs 
probability distribution is broader and shifted to heavier values, due to the bigger 
volume of parameter space allowed at higher mo and M1/2 values. Gluino mass 
distributions no longer tail off at higher masses: upper bounds would be mostly 
determined by the cut-off placed mo and Mi/2- The gluino mass and neutralino 
distributions are also much less peaked than the ones in Figs. 0D,d particularly for 
II < 0. The flatter distributions are, of course, an indication that the data aren't 
strongly constraining. Large volumes at large mo, M1/2 effectively move up the squark 
masses, as Fig. p!5| c illustrates. 

The probability distribution of Qoivih'^ is shown in Fig. [1^ for each sign of /i, when 
we drop the relic dark matter contribution. In the figure, we see that huge values 
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can result: in fact the mean value for the n > sample tails off at QoRih"^ 
much larger than the WMAP3 value of 0.1. 

This illustrates the fact that regions of param- 
eter space which fit the dark matter data in mSUGRA 
need some special annihilation mechanism that is 
not typical of the whole space. This and similar ar- 
guments have led several authors to consider non- 

where the relic density 



128, 



prior 


fiat small log 


P-/P+ 


0.45 0.43 0.19 



universal models 50, 51 



Table 5: Ratios of integrated 
probability for different signs of 
H with no dark matter con- 

might be less fine-tuned. The distributions are highly 

skewed, having tiny tails up to Qh"^ ~ 1000. The relative normalisation of /i > to 
< is calculated in the same manner as in the previous section for different priors, 
and displayed in Table ||. /i < is hardly disfavoured in R-parity violating mSUGRA 
(where we can neglect the dark matter constraint), where P_/P+ = 0.19 — 0.45. 



5. Conclusions 

We have performed global fits to mSUGRA using indirect data and state-of-the-art 
predictions of the observables. The MCMC technique was successfully employed 
despite initial non- convergence of the n < chains. Bridge sampling was used to 
normalise two isolated maxima that had not been traversed by any chain for n < 0. 
We found that /i < is somewhat disfavoured in comparison to yU > but not by huge 
margins. The rest of the fit prefers rather heavy SUSY and so the SUSY contribution 
to {g — 2)^ is small whichever the sign of /i. /i < is only disfavoured marginally, the 
ratio of integrated probability densities being P_/P+ = 0.07 — 0.16 depending upon 
the prior. We see from Fig. |l^ that without the dark matter constraint, QoMh"^ is 
predicted up to values of around 128. This corresponds to a xh of around 3x10®, 
much larger than is likely from the other observables such as M]y or sin^ (eff ) . The 
fits are therefore completely dominated by the dark matter relic density constraint 
and volume effects. Expectations that the SUSY scale will be light because of a 
preference from weak observables turn out not to be true in the global fits. If the 
dark matter constraint is dropped, as would be the case for R-parity violation, fi < 
is hardly disfavoured at all, P^/P^ = 0.19 — 0.45. fi < Ois much less disfavoured than 
many seem to assume. Many recent analyses only consider /i > on the grounds 
that /i < is strongly disfavoured by {g — 2)^. We have therefore demonstrated that 
this is not the case when one considers the entirety of the data and that /i < should 
still be considered in mSUGRA analyses. 

It could be argued that the fiat measure used here in mo, M1/2, ^0 and tan/? 
can be improved upon. For instance, tan is really a derived quantity and is related 
to more fundamental Higgs potential parameters, which could be considered more 
natural to have a fiat measure upon. There is also the issue of fine-tuning, recently 
illuminated in Ref. we could disfavour regions of parameter space that are highly 
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fine-tuned, for instance |0. Clianges sucfi as tliese in tlie prior could potentially 
change the results of the fits and we intend to investigate them in a future publication. 

Clearly, more data is required to decrease the dependence of results on the prior. 
The most helpful data is likely to be that from colliders. The MCMC fitting tech- 
nique has been used in an ATLAS study examining how cross-section and kinematic 



endpoint information constrains mSUGRA and non- universal models [53]. At this 



moment, without data from colliders, we are forced to use indirect constraints for 
the observables. However, in the future it will be desirable to predict QnAih"^ given 
some SUSY collider observables If this is in contradiction with the observed 

value from cosmology, it will point to a wrong cosmological assumption, which could 
then be changed. In order to really confirm that dark matter particles have been 
produced at colliders, one requires compatibility with direct dark matter detection 
data. Of course one would like to drop the mSUGRA assumption and perform a 
general SUSY analysis, but for this it is likely that additional data from a future 
international linear collider would be required In any case, the techniques 

investigated in this paper should prove useful for the fits. 
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A. Markov Chain Monte Carlos 

Our Markov chain consists of a list of parameter points (x*^*^) and associated like- 
lihood densities (£(x*^*))). Here, t labels the link number in the chain. Given 
some point at the end of the Markov chain (x^*-'), the Metropolis-Hastings algo- 
rithm [^, 15^, ^ requires one to randomly pick another potential point (x) (typi- 
cally in the vicinity of x*^*)) using a proposal distribution Q(x; x*^*)). There is a large 
amount of freedom in the choice of the proposal function Q, and this freedom is usu- 
ally exploited to improve the efficiency of the sampling process. In order to ensure 
that the choice of Q does not bias the final set of samples in some way, the form of 
Q is taken into account when deciding whether to accept or reject the new point. If 
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the ratio p defined by 



P 



/:(x) g(xW;x) 
/:(x(*))Q(x; xW) 



(A.l) 



is greater than one, the new point x is appended to the chain. If p is instead less than 
one, a decision must be made to determine whether to accept or reject the proposed 
point X. The rule is that acceptance of x must occur with probability p. If accepted, 
X is added to the end of the chain. If not accepted, the point x*^*^ is copied once more 
onto the end of the chain. Whichever point makes it on to the end of the chain is 
thereafter known as x'^*"''-'^^ 

As a result of following the above steps, the sampling density of points in the 
chain becomes proportional to the density of the target distribution (such as the 
posterior probability density, or the likelihood when the prior is uniform) as the 
number of links goes to infinity, under the circumstances described in Ref. [p^] . 
The Metropolis-Hastings MCMC algorithm is typically much more efficient than a 
straightforward scan for the dimensionality of input parameter space D > 3; the 
number of required steps scales roughly linearly with D rather than as a power law. 
We take the proposal function Q to be a product of Gaussian distributions along 
each dimension k = 1,2, ... ,D centred on the location of the current point along 
that dimension, i.e. a;^*'': 



where Ik denotes the width of the distribution along direction k. For the case where 
we include the dark matter relic density in the calculation, we choose Imo = 100 GeV, 
I Ml/2 ~ ^0 GeV, Iao = 400 GeV and /tan/3 = 3. For the Standard Model inputs, we 
choose Ik = 8(jfc/20. We discuss why these particular values were chosen in the next 
section. 

In order to start the chain we follow the following procedure, which finds a point 
at random in parameter space that is not a terrible fit to the data. We pick some 
y'^*') at random in the mSUGRA parameter space using a fiat distribution for its 
probability density function (pdf). The Markov chain for y is evolved through 2000 
steps. We then set x*^°^ = y^^ooo)^ continuing the Markov chain in x and discarding 
the "burn-in" chain y. A reasonable-fit point is typically found long before 2000 
iterations of the Markov chain. We must make sure that we perform enough iterations 
after this point that the chain traverses the remaining viable parameter space. We 
will provide a convergence test to this effect. 

A.l Efficiency 

The efficiency of a chain can be defined as "the number of links whose coordinates 
differ from those of their predecessor in the chain" divided by "the total number of 



D 



1 




(A.2) 
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points in the chain". There will always be a tension between efficiency and conver- 
gence^ in chains: if the Ik are set to be too small, efficiency will increase but the 
chain will take too long to achieve convergence whereas if they are too large, the 
efficiency will be so small that the sampling will contain large statistical fluctua- 
tions. In practice, the bulk of our posterior probability density is contained in a very 
thin hyper-surface in the 8-dimensional input parameter space It is thin because 
^DKih^ varies very rapidly over mSUGRA space compared to the high accuracy of 
the empirical constraint. With the Ik listed above, we found that efficiencies were at 
the per-mille level, too small to achieve a statistically stable result in a reasonable 
amount of CPU time. To achieve a significantly larger efficiency, we had to reduce 4 
to such a level that we lost convergence because the chains had not traversed the vi- 
able parameter space. In order to counter this, we expanded the errors on VLoMh^ to 
±0.02 while calculating the likelihood density in the chain. This artificially thickens 
the "surface" containing the bulk of the posterior probability density, increasing the 
efficiency to much more reasonable values of around 5-7%. In order to correct for this 
artificial thickening and re-impose the required constraint of Eq. |2.7| , it was therefore 
necessary to re- weight each link of each chain at the end of the sampling. Each link is 
re- weighted by the ratio of the proper likelihood density £p to the likelihood density 
with inflated errors £/: 



in order to impose the correct penalty on the links. We ignore additional constants 
that are independent of QjjMh'^ in this expression since the overall normalisation of 
the likelihood density is here undetermined. The re-weighting procedure necessarily 
degrades the statistical spread of the results, however we find that the increase in 
efficiency more than compensates for this effect. Below, we re- weight different vari- 
ables in order to investigate various features in the results, but the method remains 
analogous to the one described here. 

A. 2 Bridge Sampling 

One of the numbers we will require from our MCMC samples is the ratio of integrated 
posterior probabilities of yU > (-P+) and fi < (-P-). This ratio will tell us the extent 
to which /i < is disfavoured over fi > 0. Assuming a fiat prior in the variables 
of the model, the posterior probability is equal to the integrated likelihood divided 
by a factor which does not depend upon model hypotheses or parameters. Thus 
P_/P+ = / (ix£_(x)/ f (ix£+(x), where £+^_(x) is the likelihood density of /i > 
(< 0) mSUGRA respectively at parameter point (x). One way to estimate this ratio 
would be to include the sign of free parameter in the Metropolis-Hastings 

^For a discussion of convergence, see appendix H. 




(A.3) 
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procedure, to be chosen randomly in a proposal point. This algorithm leads to large 
inefficiencies because the /i > and /i < likelihood surfaces have a limited overlap, 
meaning that too many proposals for an opposite sign of /i will be rejected. Also, 
the procedure would likely provide large statistical fluctuations for the disfavoured 
sample, which we expect to be the /i < one. Even though it is disfavoured, we 
should like to investigate its properties. 

A simple way one might hope to evaluate the ratio P^/P^ is 



Pi 



1 



1 ^ 
N ^ 



(^)^ 



(A.4) 



where E_ denotes the expectation with respect to the < likelihood distribution 
[0. denotes the number of MCMC steps. Unfortunately, a simple importance 
sampling estimate of this kind does not work if there are any valid points (£ ^ 0) 
for one sign of /i that are invalid (£ = 0) for the opposite sign of /i. In our case 
there are plenty of these dangerous pairings, as sparticle mass or tachyonic bounds 
move around in parameter space depending upon the sign of [i. To get around this 



problem, we use a solution known as bridge sampling [0 with a "geometric bridge" . 
This allows us to generate a (biased) estimator r for the ratio P_/P+ as long as there 
is some viable region of /i > parameter space that is also viable for < 0. The 
estimator for the ratio is constructed as follows: 
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where (x)^^_ are the parameter points of the links in the /i > or /x < chains 
respectively. Here, we have assumed an equal number of links in each chain. In 
summary, to calculate r, we must run two chains, one for positive /i and one for 
negative n, and for every link, record the likelihood one would have obtained for 
identical input parameters except for the opposite sign of /i. 



B. Convergence and Normalisation 

In order to evaluate the convergence of the MCMC chains, we always run 9 
independent chains with different random starting points. By comparing the simi- 
larity of the resulting sampling densities of input parameters in the chains, one can 



construct a measure of convergence R. R is an upper bound on the reduction 
in variance of parameters that would result from running the chains for an infinite 
number of steps. The precise implementation is listed in Ref. j^. Values close to 1 
indicate convergence of the chains. 
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Figure 18: The two types of negative ^ samples: (a) the "/i < 0, 2 chains" samples (later 
shown to be co-annihilation samples), and (b) the "/x < 0, 7 chains" samples (later shown to 
be non co-annihilation samples, dominated by resonant Higgs annihilation regions). The 
posterior probability is indicated by the bar on the right hand side. The inner (outer) 
contours show the boundary of a 68% (95%) confidence region. 



H>0 
H<0 

|i<0, 7 chains 
|i<0, 2 chains 



We run 9 chains of 500 000 points for 
fj, > mSUGRA and for the /i < dark 
side of mSUGRA. The fi > curve in 
Fig. |l^ shows good convergence is achieved 
by 500 000 MCMC steps. However, the 
/i < curve shows a problem: conver- 
gence is never achieved. This is a seri- 
ous difficulty as one could not draw any 
quantitative statistical inferences from the 
non-converged chains. Further inspection 
of the /i < results shows that two of the 
H < chains are in a completely differ- 
ent part of parameter space than the other 
seven. This indicates isolated maxima of 
likelihood density which the MCMC has 
not been able to jump between in the fi- 
nite number of MCMC steps attempted^. There is no balance between the two 
isolated maxima in the sample. Isolating the two anomalous /i < chains and cal- 
culating R between just them, we obtain the "yU < 0, 2 chains" curve, which closely 
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step/1000 

Figure 17: Convergence statistics for the 
MCMCs. 



''A proposal distribution with longer tails, such as an n— dimensional Cauchy distribution, would 
have more chance of making such a jump. 
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approaches 1 by 500 000 MCMC steps. Thus within this isolated maximum, conver- 
gence is achieved. The same can be said of the other "/x < 0, 7 chains" samples: they 
also converge amongst themselves. Thus, the shapes of each isolated maximum are 
well determined, but the relative normalisation of the two different types of negative 
/i samples is not. 

In order to illustrate the two max- 
ima, we marginalise the two types of 
negative fi samples onto the mo — Aq 
plane in Figs. |l8|(a) and 0(b). The max- 
ima are isolated in this plane (as well as 
in some other 2-parameter planes). The 
two regions are completely separated. Theirg, 
shape is primarily determined by regions g" 
which efficiently deplete the relic den- 
sity of neutralinos which, in mSUGRA, 
is often higher than the WMAP3 con- 
straint. We investigate the 2-chain sam- 
In the figure, there are 
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Figure 19: Probability density in the light- 
est stop-lightest neutralino mass plane for the 
2-chain sample. The posterior probability is 
indicated by the bar on the right hand side. 
The inner (outer) contours show the bound- 
ary of a 68% (95%) confidence region. 



two good-fit regions: where the stau co- 
annihilates ^l|] (at moderate values of 
Ao, higher values of in the figure) 
with the LSP fiXi t7 and where the 
lightest stop co-annihilates (at Aq < —3 
TeV) with the LSP ixi tg in the 
early universe, where m^o ^ mj^, the lower strip in the figure. There was no 



significant stop co-annihilation ^ ^ region for /i > 0. On the other hand, 
Fig. |l^(b) dominantly consists of resonant Higgs annihilation regions ^ 0], 
where XiXi —>■ h,A^ ^ bb/r^r' and the focus point region where the LSP has a 
significant higgsino component and XiXi 



ZZ, WW, tt p, g|, g7 



We need a method to determine the relative normalisation of the 2-chain co- 
annihilation sample and the 7-chain resonant higgs annihilation sample. Equivalently 
we need a method to determine the ratio of the posterior probability Pi of the 2-chain 
co-annihilation sample and the posterior probability Pl of the 7-chain resonant higgs 
annihilation sample. We use P to denote the fact that the posterior probabilities are 
un-normalised. Fortunately, Eq. |A.5| provides us with a solution: we first determine 
the normalisation of the /i > sample with respect to each separate fi < sample, 
i.e. Pi/P+ and P_/P+- Since these quantities individually have good convergence 
properties, their ratio is also well determined: 



pf 

Pi 



pf 



pt_ 
pI 



0.097^0.063 = 1.53. 



(B.l) 
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Normalising the probabilities as aP^ = Pi, bPt = Pf, we fix a and b by imposing 
P^+P_ = 1 and Eq. [B.l| . The relative posterior probabilities ratios are re-calculated 
whenever alternative priors are investigated. In section |^ where we present the total 
fi < sample results, we present posterior probability densities with the correct 
normalisation, according to this prescription. The ratio of the probability of /i < 
to /i > is then determined simply by: 

p pi pf 

It is worth noting that, had we been unlucky, we might have obtained only chains 
like those in the 7-chain sample. In that case we would have carried on with the 
analysis without realising about the different 2-chain sample, therefore any results 
achieved would have been incomplete. An obvious question is: were there any other 
local maxima that we have missed by not running enough chains? Unfortunately, 
any fitting procedure is susceptible to this caveat and there is no satisfactory answer. 
Finding a high but very narrow global maximum is an unsolved problem in any 
number of dimensions. 
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